--- title: DEMS / Topography tools keywords: fastai sidebar: home_sidebar summary: "Provides a LOLA_TOPO class that manages `dem`, `slope`, and `aspect` maps." description: "Provides a LOLA_TOPO class that manages `dem`, `slope`, and `aspect` maps." nb_path: "01_dems.ipynb" ---
{% raw %}
{% endraw %} {% raw %}
{% endraw %} {% raw %}

class LOLA_TOPO[source]

LOLA_TOPO(dataset='divgdr', lat_limit=None)

Data reader class for LOLA Topography data.

It accesses preproduced data for 128 ppd for DEM, slope, and aspect,
located on the luna4 disk.

It uses virtual dask.Arrays so that virtually no memory is consumed until you
resolve a chain of operations with the `.compute()` call.

Attributes
----------
dem: xarray.DataArray
    LOLA DEM 128 ppd
slope: xarray.DataArray
    Slope in percent, derived from `dem` using `gdaldem` command line tool.
aspect: xarray.DataArray
    Aspect in degrees, derived from `dem` using `gdaldem` command line tool.
{% endraw %} {% raw %}
{% endraw %} {% raw %}
topo = LOLA_TOPO()
{% endraw %} {% raw %}
topo.dem
<xarray.DataArray 'dem' (lat: 23040, lon: 46080)>
dask.array<mul, shape=(23040, 46080), dtype=float64, chunksize=(2048, 4096), chunktype=numpy.ndarray>
Coordinates:
    y        (lat) float64 2.729e+06 2.729e+06 ... -2.729e+06 -2.729e+06
    x        (lon) float64 -5.458e+06 -5.458e+06 ... 5.458e+06 5.458e+06
  * lat      (lat) float64 90.0 89.99 89.98 89.98 ... -89.98 -89.98 -89.99 -90.0
  * lon      (lon) float64 0.0 0.007812 0.01562 0.02344 ... 360.0 360.0 360.0
Attributes: (12/19)
    transform:                  (236.901, 0.0, -5458199.04, 0.0, -236.901, 27...
    crs:                        +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=...
    res:                        (236.901, 236.901)
    is_tiled:                   1
    nodatavals:                 (-32768.0,)
    scales:                     (0.5,)
    ...                         ...
    PRODUCT_ID:                 "LDEM_128_TOPO_IMG"
    START_TIME:                 2009-07-13T17:33:17.246
    STOP_TIME:                  2011-11-14T03:44:53
    TARGET_NAME:                MOON
    long_name:                  Elevation
    units:                      m
{% endraw %} {% raw %}
topo.convert_to_lon180()
{% endraw %} {% raw %}
topo.dem
<xarray.DataArray 'dem' (lat: 23040, lon: 46080)>
dask.array<getitem, shape=(23040, 46080), dtype=float64, chunksize=(2048, 4096), chunktype=numpy.ndarray>
Coordinates:
    y        (lat) float64 2.729e+06 2.729e+06 ... -2.729e+06 -2.729e+06
    x        (lon) float64 118.5 355.4 592.3 829.2 ... -592.3 -355.4 -118.5
  * lat      (lat) float64 90.0 89.99 89.98 89.98 ... -89.98 -89.98 -89.99 -90.0
  * lon      (lon) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
Attributes: (12/19)
    transform:                  (236.901, 0.0, -5458199.04, 0.0, -236.901, 27...
    crs:                        +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=...
    res:                        (236.901, 236.901)
    is_tiled:                   1
    nodatavals:                 (-32768.0,)
    scales:                     (0.5,)
    ...                         ...
    PRODUCT_ID:                 "LDEM_128_TOPO_IMG"
    START_TIME:                 2009-07-13T17:33:17.246
    STOP_TIME:                  2011-11-14T03:44:53
    TARGET_NAME:                MOON
    long_name:                  Elevation
    units:                      m
{% endraw %} {% raw %}
topo.convert_to_lon360()
{% endraw %} {% raw %}
topo = LOLA_TOPO(lat_limit=80)
{% endraw %}

{% include note.html content='The keys of the fnames dict become attributes in the class, storing the opened xarrays.' %} The xarray.DataArray object is a very rich object for n-dimensional datacubes with physical dimensions.

{% raw %}
topo.dem
<xarray.DataArray 'dem' (lat: 20480, lon: 46080)>
dask.array<getitem, shape=(20480, 46080), dtype=float64, chunksize=(2048, 4096), chunktype=numpy.ndarray>
Coordinates:
    y        (lat) float64 2.426e+06 2.426e+06 ... -2.426e+06 -2.426e+06
    x        (lon) float64 -5.458e+06 -5.458e+06 ... 5.458e+06 5.458e+06
  * lat      (lat) float64 80.0 79.99 79.98 79.98 ... -79.98 -79.98 -79.99 -80.0
  * lon      (lon) float64 0.0 0.007812 0.01562 0.02344 ... 360.0 360.0 360.0
Attributes: (12/19)
    transform:                  (236.901, 0.0, -5458199.04, 0.0, -236.901, 27...
    crs:                        +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=...
    res:                        (236.901, 236.901)
    is_tiled:                   1
    nodatavals:                 (-32768.0,)
    scales:                     (0.5,)
    ...                         ...
    PRODUCT_ID:                 "LDEM_128_TOPO_IMG"
    START_TIME:                 2009-07-13T17:33:17.246
    STOP_TIME:                  2011-11-14T03:44:53
    TARGET_NAME:                MOON
    long_name:                  Elevation
    units:                      m
{% endraw %}

Slicing

Constraining latitude

Sometimes we need to cut off some latitudes because some products have been produced only up to a certain latitude.

This is necessary if one needs the exact same pixels as another map product, for pixel-based slicing through multiple data products. Otherwise, xarray offers also direct index-label-value-based access, i.e. lat/lon coordinates.

{% raw %}

LOLA_TOPO.slice_lat[source]

LOLA_TOPO.slice_lat(data, lat)

Return the map `data` constrained to lat <= `lat`.

Parameters
----------
data: {'dem', 'slope','aspect'}
    String that choses which data product should be constrained.
lat: int, float
    Limiting latitude value.
{% endraw %} {% raw %}
topo.slice_lat("slope", 70)
<xarray.DataArray 'slope' (lat: 17920, lon: 46080)>
dask.array<getitem, shape=(17920, 46080), dtype=float64, chunksize=(2048, 4096), chunktype=numpy.ndarray>
Coordinates:
    y        (lat) float64 2.123e+06 2.122e+06 ... -2.122e+06 -2.123e+06
    x        (lon) float64 -5.458e+06 -5.458e+06 ... 5.458e+06 5.458e+06
  * lat      (lat) float64 70.0 69.99 69.98 69.98 ... -69.98 -69.98 -69.99 -70.0
  * lon      (lon) float64 0.0 0.007812 0.01562 0.02344 ... 360.0 360.0 360.0
Attributes: (12/19)
    transform:                  (236.901, 0.0, -5458199.04, 0.0, -236.901, 27...
    crs:                        +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=...
    res:                        (236.901, 236.901)
    is_tiled:                   1
    nodatavals:                 (-32768.0,)
    scales:                     (0.01,)
    ...                         ...
    PRODUCT_ID:                 "LDEM_128_SLOPE_IMG"
    START_TIME:                 2009-07-13T17:33:17.246
    STOP_TIME:                  2011-11-14T03:44:53
    TARGET_NAME:                MOON
    long_name:                  Slope
    units:                      deg
{% endraw %}

Rectangular slicing

Note, reading with chunks implies lazy execution, meaning to get actual values, .compute() needs to be called.

{% raw %}
topo.get_slice("slope", 20, 131, dlat=1, dlon=1)
<xarray.DataArray 'slope' (lat: 128, lon: 129)>
dask.array<getitem, shape=(128, 129), dtype=float64, chunksize=(128, 129), chunktype=numpy.ndarray>
Coordinates:
    y        (lat) float64 6.367e+05 6.364e+05 6.362e+05 ... 6.068e+05 6.066e+05
    x        (lon) float64 -1.486e+06 -1.485e+06 ... -1.456e+06 -1.455e+06
  * lat      (lat) float64 21.0 20.99 20.98 20.97 ... 20.03 20.02 20.01 20.0
  * lon      (lon) float64 131.0 131.0 131.0 131.0 ... 132.0 132.0 132.0 132.0
Attributes: (12/19)
    transform:                  (236.901, 0.0, -5458199.04, 0.0, -236.901, 27...
    crs:                        +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=180 +x_0=...
    res:                        (236.901, 236.901)
    is_tiled:                   1
    nodatavals:                 (-32768.0,)
    scales:                     (0.01,)
    ...                         ...
    PRODUCT_ID:                 "LDEM_128_SLOPE_IMG"
    START_TIME:                 2009-07-13T17:33:17.246
    STOP_TIME:                  2011-11-14T03:44:53
    TARGET_NAME:                MOON
    long_name:                  Slope
    units:                      deg
{% endraw %} {% raw %}
azi_slice = topo.get_slice("aspect", 20, 131, dlat=1, dlon=1)
azi_slice.compute()
<xarray.DataArray 'aspect' (lat: 128, lon: 129)>
array([[139. , 139.5, 135.9, ...,  27.6,  16.6,   8.1],
       [143.4, 141.7, 136.1, ...,  37. ,  20.6,  14.5],
       [144.7, 141.1, 135.4, ...,  39.2,  28.4,  25.2],
       ...,
       [315.6, 333.2, 307.8, ..., 134.2, 132.3, 131.5],
       [299.7, 306.5, 312.6, ..., 137.5, 134.2, 131.9],
       [321.4, 322.4, 316.1, ..., 141.6, 140.1, 135.9]])
Coordinates:
    y        (lat) float64 6.367e+05 6.364e+05 6.362e+05 ... 6.068e+05 6.066e+05
    x        (lon) float64 -1.486e+06 -1.485e+06 ... -1.456e+06 -1.455e+06
  * lat      (lat) float64 21.0 20.99 20.98 20.97 ... 20.03 20.02 20.01 20.0
  * lon      (lon) float64 131.0 131.0 131.0 131.0 ... 132.0 132.0 132.0 132.0
Attributes:
    long_name:  Aspect
    units:      deg(north)
{% endraw %} {% raw %}
assert isinstance(topo.dem_fpath, Path)
assert isinstance(topo.slope_fpath, Path)
assert isinstance(topo.aspect_fpath, Path)
assert isinstance(topo.dem, xr.DataArray)
assert isinstance(topo.slope, xr.DataArray)
assert isinstance(topo.aspect, xr.DataArray)
{% endraw %} {% raw %}
assert topo.dem.ndim == 2
assert topo.slope.ndim == 2
assert topo.aspect.ndim == 2
{% endraw %}

Longitudes

Many map-based applications still use the 180 degree longitude system, but some newer ones use Lon360.

There are converter methods that switch back and forth between these 2 coordinate systems.

{% raw %}

LOLA_TOPO.convert_to_lon180[source]

LOLA_TOPO.convert_to_lon180()

Switch all three maps to -180/180 longitude system.
{% endraw %} {% raw %}
topo.convert_to_lon180()
{% endraw %} {% raw %}
topo.dem.lon
<xarray.DataArray 'lon' (lon: 46080)>
array([-180.      , -179.992188, -179.984375, ...,  179.976562,  179.984375,
        179.992188])
Coordinates:
    x        (lon) float64 118.5 355.4 592.3 829.2 ... -592.3 -355.4 -118.5
  * lon      (lon) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
{% endraw %} {% raw %}

LOLA_TOPO.convert_to_lon360[source]

LOLA_TOPO.convert_to_lon360()

Switch all three maps to 360 longitude system.
{% endraw %} {% raw %}
topo.convert_to_lon360()
{% endraw %} {% raw %}
topo.dem.lon
<xarray.DataArray 'lon' (lon: 46080)>
array([0.000000e+00, 7.812500e-03, 1.562500e-02, ..., 3.599766e+02,
       3.599844e+02, 3.599922e+02])
Coordinates:
    x        (lon) float64 -5.458e+06 -5.458e+06 ... 5.458e+06 5.458e+06
  * lon      (lon) float64 0.0 0.007812 0.01562 0.02344 ... 360.0 360.0 360.0
{% endraw %}

Plotting methods

{% raw %}
demplot = topo.plot_dem(20, 131, dlon=2)
{% endraw %} {% raw %}
demplot
{% endraw %} {% raw %}
slopeplot = topo.plot_slope(20, 131, dlon=2)
{% endraw %} {% raw %}
(demplot + slopeplot).cols(1)
{% endraw %} {% raw %}
topo.plot_aspect(20, 131, dlon=2)
{% endraw %} {% raw %}
lat = 20.572
lon = 131.301
{% endraw %} {% raw %}
topo.get_elev_by_coord(lat, lon)
263.25
{% endraw %} {% raw %}
topo.get_az_by_coord(lat, lon)
197.1
{% endraw %}